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We describe how to implement the time-dependent variational principle for matrix product states in the 
thermodynamic limit for nonuniform lattice systems. This is achieved by confining the nonuniformity to 
a (dynamically growable) finite region with fixed boundary conditions. The suppression of unphysical 
quasiparticle reflections from the boundary of the nonuniform region is also discussed. Using this algo- 
rithm we study the dynamics of localized excitations in infinite systems, which we illustrate in the case of 
the spin-1 anti-ferromagnetic Heisenberg model and the 0^ model. 



Douglas Adams (nearly) put it best: "[Hilbert] space is 
big. ... You just won't believe how vastly hugely mind- 
bogglingly big it is. I mean, you may think it's a long way 
down the road to the chemist, but that's just peanuts com- 
pared to [Hilbert] space." Given said space's exponential 
growth with the size of a many-particle system, it is a little 
astounding that general techniques exist to allow efficient 
numerical calculations in a wide range of physically inter- 
esting cases. This is possible because physically relevant 
states have limited entanglement [1-3]. This observation 
may be exploited to obtain an efficient parametrization of 
this physical corner of Hilbert space. 

The class of matrix product states (MPS) [4J represent, 
in one dimension, a good parametrization of the physi- 
cal corner. This is amply demonstrated by the unparal- 
leled success of the density matrix renormalization group 
(DMRG) [5 1, which can be viewed as a variational method 
applied to MPS [6|. The MPS class has served as the ba- 
sis for many exciting generalisations, including the study 
of non-equilibrium dynamics fT| and higher-dimensional 
systems [8|. More recently, Haegeman et al have imple- 
mented the time-dependent variational principle (TDVP 
— see boxout) for MPS [9 |, providing a locally optimal (in 
time) framework for finding ground states, simulating dy- 
namics, and studying excitations of one-dimensional lattice 
systems. 

The simulation of infinite quantum spin systems has 
been mostly confined to the translation invariant setting 
(usually by restricting states to translation invariant subsets 
of MPS). However, the ability to explore locally nonuni- 
form states on an infinite lattice is particularly attractive 
for studying the dynamics, e.g. scattering, of localized ex- 
citations in large systems. For example, this would provide 
a realistic setting in which to study quantum field excita- 
tions. There has been some prior work in this direction, 
building on the light-cone results of EKTOl, where the dy- 
namics of a local disturbance is (partially) studied in the 
Heisenberg picture. These approaches can become expen- 
sive for systems with large local spin dimensions (such 



as those appearing in lattice field theory). Another di- 
rection, suggested in [11 J, is to work completely in the 
Schrodinger picture with infinite uniform MPS and to add 
a finite nonuniform region. 

In this Letter we explore the locally optimal implemen- 
tation of the TDVP for uniform MPS with a dynamically 
growable nonuniform segment. We derive the equations 
of motion for the variational parameters using a partic- 
ular choice of gauge-fixing which allows us to integrate 
the variational dynamics with a complexity that scales as 



The time-dependent variational principle 

A variational mani- 
fold M is depicted as 
embedded in a Hilbert 
space H. Beginning with 
a state |^[a(t = 0)]) in 
M, where a(t) are the 
variational parameters, 
we wish to compute the 
time-evolution accord- 
ing to the Schrodinger 
equation ^ 




d \na{t)\) 



The exact evolution generally leads out of M. Equiv- 
alently, the infinitesimal time step —\H\^[a{t)]) (the blue 
arrow) need not lie within the tangent plane T to y\4 at 
point |^[a(t)]) (the green dashed line). The best approxi- 
mation to the exact evolution, whilst remaining in re- 
quires a tangent vector |$) G T (the red arrow) that best 
approximates —iH\^[a{t)]), which is found by projecting 
—\H \^f[a(t)\) onto T. In other words, |$) must minimize 
||iH|*[a(t)]) + |$)||2. 

This is equivalent to finding optimal equations of mo- 
tion for a, since writing |$) = (where := 
d/da^ |^[a(t)])) and taking the derivative of the above mag- 
nitude with respect to a"^ results in the flow equation 

where g^^ is the inverse of gj]^ — {dj^\d]^^), which is the 
pullback metric on T. Here, we assume that |^[a(t)]) is a 
holomorphic function of a(t), although this is not necessary. 
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d\t\D^N, where is the length of the generic piece (the 
number of sites), |t| is the desired integration time, d is the 
local spin dimension, and D is the bond dimension. Even 
though the ends of the generic region can move, there may 
be some backscattering due to boundary effects; we de- 
scribe how to compensate for these with the addition of 
an optical potential term. This method is illustrated in the 
case of local excitations of the spin-1 anti-ferromagnetic 
Heisenberg model and for particles in (j)^ theory. 

We assume throughout that our Hamiltonian H con- 
tains only nearest-neighbour terms. It is decomposed as 
H = i/^^i + where i/^^^ = E^^-oo K^n+i with 
Kn^i = htm^,. yn^m, and H'^^ = EnJi' <^n+i 
with [1, A^] representing a compact contiguous region of 
the lattice and h^n%+i = forn < 1, n > A^, allowing us 

to also write if = E„ K,n+i = E„ K^+i + h'°%+,] . 
We consider two cases in particular: Firstly, a non-trivial 
h^^^ leads to a locally nonuniform ground state, which can 
be found using imaginary time evolution via our algorithm. 
Secondly, given a purely uniform Hamiltonian (h^^^ = 0) 
and an initial state that differs only locally (in a region 
[1, A^]) from an eigenstate of H^^\ our algorithm can be 
used to simulate the resulting locally non-trivial dynamics. 

To capture a locally nonuniform state using MPS, we de- 
fine a class of "sandwich" states (sMPS), based on uniform 
MPS, using two dx D x D tensors Al and Ar describing 
the (asymptotic) state either side of the nonuniform region 
[1, A^], which is parametrized by A^ further tensors. An 
sMPS state can be written as 
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where |s) = |. . . 5i . . . . . .) and A^^ G Md(C) (where 

X = L,R, [1, A^]). Taking Al = Ai = • • • = An = 
Ar gives a completely uniform state. The vectors Vl/r 
are, as with uniform MPS [9|, generically irrelevant to the 
TDVP algorithm and are not further specified. In principle, 
the dimensions of A^ are subject only to the constraints 
of the matrix product, which can become important when 
maximizing numerical efficiency. However, for reasons of 
notational simplicity, we assume uniform dimensions here. 

Al/r represent the left and right asymptotic states: The 
reduced density matrix P[n,m](^L, Ar, ^i . ^r) of a piece 
of the lattice in the left or right region n,m < lorn,m > 
A^ tends to that of the uniform MPS state P[n,m](AL/R) as 
the distance from the nonuniform region increases. Since 
Al /r represent infinite "bulk" regions of the lattice, they 
should not be affected by nonuniformities in the [1, A^] re- 
gion. Further, if the left and right asymptotic states are 
eigenstates of H^^\ they are left completely unchanged 
by time evolution. Assuming this, we restrict the varia- 
tional parameters to the tensors Ai . . . A^ and treat Al/r 
as boundary conditions, resulting in NdD'^ parameters for 
the sMPS variational manifold M (1%. Al/r can be ob- 
tained for the ground state of H^^' using the existing TDVP 



algorithm for uniform MPS. To accurately capture states 
with a nonuniform region [1, A^] in this way, A^ should be 
sufficiently large so that the asymptotic states are already 
reached at the left and right boundaries with the bulk. 

The tensor network formed by the matrices A can be 
visualized as 

(WkW>4<^^4^4<^Wk;k:>---- 

with the nonuniform region in the centre and the physical 
indices pointing upwards. 

A state so represented is invariant under gauge transfor- 
mations 

Al QlAIqI^ 

Ar ^ 9rA19r^ 

with gx G Md{C), and where Qq = Ql and = Qr. 
Since we allow only jv to vary, the tangent plane T 
to M contains (A^ — 1)D'^ — 1 gauge degrees of free- 
dom corresponding to ^i...Ar-i, excluding the trivial gauge 
transformation gi,,,N-i ^ I- The dimension of A4 is thus 
reduced from A^dZ)2 todim(7W) = {N{d-1) + 1)D^+1. 
Restricting to normalized states further reduces the dimen- 
sion by 1. This redundancy of parameters will be used to 
simplify the TDVP algorithm. 

As with the uniform MPS TDVP algorithm, we re- 
quire Al/r such that the transfer matrices El/r = 

J2t^L/R ® ^1/R h^^^ 1 as a unique eigenvalue with 
largest absolute value. We denote the corresponding left 
and right eigenvectors as {Il/r\ , |^L/i?), with correspond- 
ing matrix representations Il/r^t^l/r- We further define 

ln>i Yfs^n^n-iAl with /^<i = II and rn<N 

with rn>N = T whcrC A^yN = 

Ar and A^^i = Al^ The expectation value of an operator 
On acting on a single site n anywhere on the chain is then 

{on) = (^n-i| {s\on\t) A^ A^ [r^), wMch can be 
extended easily for more general operators. 

To implement the TDVP (see boxout), \\\H \^[A{t)]) ^ 
1$) IP must be minimized with respect to the tan- 
gent vector 1$) G T, which can be written as 

MB]) = T.nS=lBn.^ l^n,.*[^]) (with \dnMA]) = 

d/dAn^i |^^[^])). Expanding leaves terms {^[B]\^[B]) 
and ($[5] |i7 + h.c, where the remaining term 

is a constant that can be ignored. The metric ($[5] 
is, at first glance, very complicated since it couples the 
tensors for different lattice sites, precluding a split- 
ting of the problem into A^ separate parts (one for each 
Bn). Fortunately, as in [9 1, it is possible to use the gauge 
freedom in the tangent plane to impose gauge-fixing con- 
ditions (GFC) Yfs K ® K K) = 0, Vn Gl2, A^], such 
that all mixing terms are eliminated and = 
^n=i i^n-i IZlf ® ^1 l^n)- Wc therefore parametrize 
the B^n=2...N t>y defining the {d — 1)D x dD matrix 
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to contain an orthonormal basis for the null space 

of [Rn\(a,sy,l3 = k^^^^^la,/? SO that ^{Xn) = 

ln-i^V^''^~^^'^ always satisfies the GFC and finally 

The GFC and corresponding parametrization for B2...N 
are the same as for the MPS case, albeit using the above 
extended definitions for 1^ and r^. As in [9|, it is also 
straightforward to prove that imposing the GFC uses up 
exactly the gauge degrees of freedom, the key difference 
being that, in this case, there is not enough gauge free- 
dom to apply a GFC to Bi. This does not significantly 
complicate matters, since the GFCfor B2...N already elim- 
inates all mixing terms from It does, how- 
ever, make it necessary to explicitly eliminate the norm- 
changing degree of freedom from the tangent space, be- 
cause = {lL\YfsBl®A\\ri) + 0. This 
can be done by replacing B. with B = B — (^\B\^) in 
the TDVP equations, effectively projecting out the norm- 
changing component of i7 |^^) (see [9|). 

Using the above, ($[5] also simplifies, but 

still contains terms mixing B^ and /im,m+i = /im,m+i — 
{hm^m+i) for n, n+l^m. Each B^ term con- 
tains a sum over hm,m+i extending into the right-hand 
bulk: this is understood by defining the effective hamil- 

tonian = j:Z=n+i (llZn+i E,) 

\Kn) En \K, 



(also giving 
where E^^ - 



n+l/ 



H 



' m+1/ 



n+1 



with C^'* = 
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FIG. 1 . Simulated real time evolution (using a fixed nonuniform 
region [-120, 80] with D = 64 and dt = 0.1) of the spin-1 anti- 
ferromagnetic Heisenberg model with two localized entangled 
excitations generated by applying S^_-S^j^- at m = ±15 with 
j = 5. The plots show the expectation value of Sn with the top- 
right plot showing the excitation bouncing at the right boundary. 
For the bottom-right plot, we used an optical potential — ie/in^^^+i 
rising linearly from eeo = at site 60 to ego = 1 at site 80 to sup- 
press this reflection, albeit imperfectly. The uniform ground state 
was converged up to a state tolerance 77 = 10~^^. For the time- 
integration we used a 4th order explicit Runge-Kutta algorithm. 



En,^=i {s, t\hn,n+i \u, v) A^A^^^. The uniform part of 
the sum \Kji) = \Kn^i) can be computed by exploiting 
the spectrum of Er, giving 1^^^) = (I - EnfEf \rR) 
or, equivalently, {1-Er^- |r^) (/^|) \Kr) = E^ |r^), 
which can then be solved for in the matrix represen- 
tation {0{D^)). The Bl term additionally contains an in- 
finite sum over the left-hand bulk {Kl \ which can be com- 
puted analogously to | K^) . Note that the energy difference 
due to the nonuniformity is AE" = (Kx.|ro) + {In\Kr) + 

En=o (^n,n+i) " (A^ + l) {h)^^-, whcrc {h)^^- is the en- 
ergy per-site of the uniform bulk state. 

We now have the ingredients needed to implement the 
TDVP efficiently. To this end we define 



=2. ..AT 



E/V2 ^s,t At t -1/2t^s 



s,t 



s,t 
d 



and 



s 
d 



that mB]\BmA]) = En=2tr[xiF.] + 

^itr iLGlnBf . With the result for ($[5] I 

the TDVP gives us — 1 + d independent matrix 
equations, which are solved when the parameter ma- 
trices = -iFr, for n G [2,7V] and Bf = -iG^ 
giving us the time evolution of the variational parameters 



so 



At 



=2. ..AT 



= Al = -iGl{A), These 



equations can be integrated numerically, for example with 
the following simple algorithm implementing the Euler 
method: 

1. Calculate F^, forn = 1 . . . TV, 5 = 1 . . . d. 

2. Take a step by setting An{t + dt) = An{t) + dtBn^ 

3. Restore a canonical form of A[i,Ar] using a gauge 
transformation and normalize the state by rescaling 
Al and Ar, 

4. Compute desired quantities, such as the energy, and 
adjust the step size dt as required. 

5. If needed, expand the nonuniform region to the left 
and/or right. 

Normalization is necessary because the norm is only pre- 
served to first order in dt. Maintaining a canonical form 
(e.g. To. ..AT = I) can simplify some parts of the TDVP 
calculations and improve the conditioning of the matrices 



4 



O 

X 



12001 
1000 1 
800 
600 1 
400 1 
200 1 

o' 



I 



■ 400 800 1200 
time(xO.lfes) 



-300 -200 -100 100 200 300 
site 



3.25 
3.00 
2.75 
2.50 
2.25 
2.00 
1.75 
1.50 
1.25 
1.00 



(a) 



(b) 



(c) 



FIG. 2. The time evolution of the block entropy S of one half 
of the lattice, split at each site n, for the same simulation as in 
Fig.[T] except that dynamic expansion of the nonuniform region 
is used and D = 16. The inset shows a cross-section at site 
(blue line). The green line shows the block entropy of the uni- 
form ground state. Note that the entanglement appears to tend to 
the asymptotic value of approximately 2. This strongly suggests 
that a hybrid method whereby uniform matrices are reintroduced 
between the two excitations as they become separated will allow 
the study of the asymptotics of entangled excitations for large 
times (this approach will be the subject of a forthcoming paper). 



involved. The last step allows for a small initial nonuni- 
form region, which can be grown if the dynamics warrant 
changing the state significantly outside of it. This is done 
by "absorbing" sites from the uniform region(s) into the 
nonuniform region, copying the and Aji matrices as 
needed. 

Whether it is necessary to grow the nonuniform region 
can be heuristically determined by comparing the contri- 
bution rjn = ^J {ln-i \ Ylt ® Bn\rn) of the individual 
Bn to the optimal tangent vector obtained above 

with the corresponding final 77uni obtained from the uniform 
MPS TDVP algorithm in finding the bulk ground state. If 
Vn > Vuni + 5 for n near the boundaries, the nonuniform 
region should be expanded until this is no longer the case. 
Note that, due to the particular (right) GFC used here, the 
errors incurred near the left boundary are different to those 
due to the right boundary. This asymmetry can be seen in 
the form of A^n+i^ which only includes the effects of the 
Hamiltonian terms to the right of n. In practice, represent- 
ing a given non-uniformity localized near a site m gener- 
ally requires the nonuniform region to extend further to the 
left than to the right of m. 

Note also that the above algorithm is not well suited to 
simulating real-time dynamics because errors due to the 
simple integration method used are cumulative. Instead, 
more sophisticated integrators such as the commonly used 
fourth-order explicit Runge-Kutta method [ 14J are prefer- 
able. The Euler method is, however, still useful for finding 
ground states because imaginary time evolution is "self- 
correcting" — it will always take you towards the ground 
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FIG. 3. Simulation (using dynamic expansion of the nonuniform 
region, with D = 24, d = 16, dt = 0.03) of particle scattering 
in real scalar field theory for lattice spacings a^X of 0.1 (a), 
0.2 (b) and 0.5 (c). Dynamical expansion of the generic region is 
indicated by the staircase pattern. The Hamiltonian if^^^ = 

En (^^n + 4^0n + + ^(0n - 0n+l)'), where 

[(/)n(t), 7rm(t)] = iSnm, IS a spacc-discretizcd version of H — 
J dx [^{7t{x, if + (V(/)(x, t)f + t)) + t)4] in 

(1 + 1) dimensions. For certain values of iiQ,\,a the ground 
state spontaneously breaks the Z2 symmetry cj) — —cj) such that 
(0) = ±00- The above plots show the evolution of two excita- 
tions created by applying e~^^^/^^ to the approximate ground 
state of H^^^ with a fixed dimensionless parameter A//i^ — 99, 
which is on the broken side of the critical point separating the 
broken and unbroken phases, ji is the physical mass related to 
the bare mass /xq via /xq = — ^M^^ where 5 11^ is determined 
by the one-loop correction, which is divergent in the continuum. 
For more details about the application of MPS to real scalar cj)^ 
theory and its critical behaviour, see 1 12| and 1 13 J . 




state, given that the starting point is not orthogonal to it. 

To demonstrate our algorithm, we use the anti- 
ferromagnetic spin-1 Heisenberg model h^^j^-^ — 

= JT.i=x,y,zS'nSl^+i' Having found a uni- 
form MPS approximation for the ground state, we use 
imaginary time evolution to solve an impurity prob- 
lem h}^^-^ = \h^™ with all other /il,^^o,n+i = 0' 
where we have chosen the 
nonuniform region to be 
[-41,33]. The plot to the 
right shows {S^) in units of 
J near site for A = — 2J 
with bond dimension D — 
and a "state tolerance" 

As a demonstration of real-time evolution, we use the 
same /i^"^ as above but with h}^^^^i — 0, Vn, and intro- 
duce two local excitations to the uniform ground state ap- 
proximation. We do this by applying the (nonunitary) op- 
erator S^_jS^^j, generating an entangled excitation, to 
two separated sites at m = i/c, with k + j < N/2. By 
calculating the expectation value of an observable such as 
for a set of sites (possibly extending into the left and 
right bulk regions) after each step, the time evolution of 
an operator can be visualized as in Fig.[T] To mitigate non- 
physical reflections that can occur when a travelling excita- 
tion meets a boundary with the uniform region, "optical po- 
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tential" terms h^^^^^i — ~i^n^n"n+i t)e locally turned 
on near to the boundaries. The parameters can be tuned 
in order to dissipate energy heading out of the nonuniform 
region. This effect can also be seen in Fig.[T] 

As a final demonstration of our approach we simulate the 
scattering of localised excitations in (j)^ theory on the lat- 
tice, depicted in Fig. [3] A sequence of results approaching 
the continuum limit is presented, demonstrating the con- 
vergence of the algorithm for small lattice spacing. 

In this paper, we have introduced an efficient means of 
simulating the dynamics of localized nonuniformities on 
spin chains in the thermodynamic limit using the time- 
dependent variational principle (TDVP) and a special class 
of matrix product states (MPS). As with the existing algo- 
rithms implementing the TDVP for MPS in other settings 
[|91, this algorithm approximates exact time evolution op- 
timally given the restrictions of the variational class. Our 
(open source) implementation [15] is available as Python 
(http://www.python.org) source code, including example 
simulation scripts. 

During completion of this work, we learned of the inde- 
pendent results of V. Zauner, M. Ganahl, H.G. Evertz, and 
T. Nishino (see arXiv: xxxx.xxxx), and H. N. Phien, G. 
Vidal, and I. P. McCulloch (see arXiv: xxxx.xxxx). 
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